The DEA water observation product maps the presense of surface water from Landsat imagery. An analysis through time can map the frequency a pixel is inundated by water, which can be used to infer the temporal and spatial statistics of flood/drought events. This notebook demonstrates how to run the DEA water observation algorithm for a given area of interest
This notebook requires the "WOfS Environment", which needs to be installed on the first run of this notebook by executing the following cell. Once done, set the environment in the top-right corner.

# Install the WOfS environment
# - Only installs the environment if required (run it at least once)
# - Select the "WOfS Environment" kernel as per the picture above
!../tools/install_wofs.sh
WOfS Environment kernel already installed
aoi_label = "Lake Tempe, Indonesia"
product = "landsat8_c2l2_sr"
longitude = (119.8242517, 120.0350519)
latitude = (-4.2013799, -3.9445384)
time = ('2020-01-01', '2020-12-31')
output_crs = "EPSG:32650"
resolution = (30, -30)
# Where to save the DEM fetched in ODC
DEM_PATH = "/home/jovyan/dems/srtm_lake_tempe.tif"
from os import environ
environ["AWS_HTTPS"] = "NO"
environ["GDAL_HTTP_PROXY"] = "easi-caching-proxy.caching-proxy:80"
print(f'Will use caching proxy at: {environ.get("GDAL_HTTP_PROXY")}')
Will use caching proxy at: easi-caching-proxy.caching-proxy:80
import sys
from pathlib import Path
import xarray as xr
import rioxarray
from wofs.virtualproduct import WOfSClassifier
try:
from dea_tools.plotting import display_map, rgb
except ModuleNotFoundError:
# Local copy of selected dea_tools
if '../tools' not in sys.path:
sys.path.append('../tools')
from datacube_utils import display_map
rgb = None # Not copied or adapted yet
# Ignore warnings in output
import warnings
from sqlalchemy.exc import SAWarning
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=SAWarning)
from os import environ
from cartopy.crs import PlateCarree
from datashader import reductions
import holoviews as hv
import hvplot.xarray
import matplotlib.pyplot as plt
from datacube import Datacube
from datacube.utils.aws import configure_s3_access
configure_s3_access(
aws_unsigned=False,
requester_pays=True,
);
dc = Datacube()
# Display available products
products_info = dc.list_products()
products_info
| name | description | license | default_crs | default_resolution | |
|---|---|---|---|---|---|
| name | |||||
| copernicus_dem_30 | copernicus_dem_30 | Copernicus 30m Digital Elevation Model (GLO-30) | None | None | None |
| landsat5_c2l2_sr | landsat5_c2l2_sr | Landsat 5 Collection 2 Level-2 Surface Reflect... | None | None | None |
| landsat5_c2l2_st | landsat5_c2l2_st | Landsat 5 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| landsat7_c2l2_sr | landsat7_c2l2_sr | Landsat 7 USGS Collection 2 Surface Reflectanc... | None | None | None |
| landsat7_c2l2_st | landsat7_c2l2_st | Landsat 7 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| landsat8_c2l1 | landsat8_c2l1 | Landsat 8 Collection 2 Level-1 (top of atmosph... | CC-BY-4.0 | None | None |
| landsat8_c2l2_sr | landsat8_c2l2_sr | Landsat 8 Collection 2 Surface Reflectance, pr... | None | None | None |
| landsat8_c2l2_st | landsat8_c2l2_st | Landsat 8 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| landsat9_c2l1 | landsat9_c2l1 | Landsat 9 Collection 2 Level-1 (top of atmosph... | CC-BY-4.0 | None | None |
| landsat9_c2l2_sr | landsat9_c2l2_sr | Landsat 9 Collection 2 Surface Reflectance, pr... | CC-BY-4.0 | None | None |
| landsat9_c2l2_st | landsat9_c2l2_st | Landsat 9 Collection 2 Level-2 UTM Surface Tem... | CC-BY-4.0 | None | None |
| lpdaac_nasadem | lpdaac_nasadem | NASADEM Merged DEM Global 1 arc second V001 | None | None | None |
| nasa_aqua_l2_oc | nasa_aqua_l2_oc | NASA MODIS-Aqua L2 Ocean Color, regridded to W... | None | EPSG:4326 | (0.01, 0.01) |
| nasa_aqua_l2_sst | nasa_aqua_l2_sst | NASA MODIS-Aqua L2 Sea Surface Temperature, re... | None | EPSG:4326 | (0.01, 0.01) |
| s2_l2a | s2_l2a | Sentinel-2a and Sentinel-2b imagery, processed... | None | None | None |
display_map(x=longitude, y=latitude)
Load the SRTM DEM data from ODC. If the area of interest is reasonably flat (like Lake Tempe) the DEM can be optional.
In this notebook, the DEM data is retrieved and saved to the same local file each time the next cells are run. This if fine if the area of interest is being changed (as per this demonstration notebook). For other use-cases, consider retrieving and saving the DEM data once only.
dem = dc.load(
product="lpdaac_nasadem",
latitude=latitude,
longitude=longitude,
output_crs="epsg:4326",
resolution=(-1/3600, 1/3600),
)
elevation = dem.elevation
options = {
'title': 'Elevation',
'width': 800,
'height': 500,
'aspect': 'equal',
'cmap': plt.cm.terrain,
'clim': (0, elevation.max().values.item()), # Limit the color range depending on the layer_name
'colorbar': True,
'tools': ['hover'],
}
plot_crs = PlateCarree()
elevation.hvplot.image(
x = 'longitude', y = 'latitude', # Dataset x,y dimension names
crs = plot_crs,
rasterize = True, # If False, data will not be reduced. This is slow to load but all data is loaded.
aggregator = reductions.mean(), # Datashader calculates the mean value for reductions (also first, min, max, las, std, mode)
precompute = True, # Datashader precomputes what it can
).opts(**options).hist(bin_range = options['clim'])
dem_path = Path(DEM_PATH)
dem_path.parent.mkdir(parents=True, exist_ok=True)
elevation.rio.to_raster(dem_path)
Load the data from ODC and rename bands as needed by the WOfS classifier.
# Ignore SAWarning in output
measurements = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
data = dc.load(
product=product,
longitude=longitude,
latitude=latitude,
time=time,
output_crs=output_crs,
resolution=resolution,
measurements=measurements,
dask_chunks={'time': 1},
)
data = data.rename({
"blue": "nbart_blue",
"green": "nbart_green",
"red": "nbart_red",
"nir": "nbart_nir",
"swir1": "nbart_swir_1",
"swir2": "nbart_swir_2",
"pixel_qa": "fmask",
})
data
<xarray.Dataset>
Dimensions: (time: 19, y: 951, x: 785)
Coordinates:
* time (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-1...
* y (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
* x (x) float64 8.371e+05 8.37e+05 ... 8.136e+05 8.136e+05
spatial_ref int32 32650
Data variables:
nbart_blue (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_green (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_red (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_nir (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_swir_1 (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
nbart_swir_2 (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
fmask (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
crs: EPSG:32650
grid_mapping: spatial_refRun the water observations classifier. This may take a few minutes.
transform = WOfSClassifier(c2_scaling=True, dsm_path=DEM_PATH)
# Compute the WOFS layer
wofl = transform.compute(data)
wofl
<xarray.Dataset>
Dimensions: (y: 951, x: 785, time: 19)
Coordinates:
* y (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
* x (x) float64 8.371e+05 8.37e+05 8.37e+05 ... 8.136e+05 8.136e+05
* time (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-12...
spatial_ref int32 32650
Data variables:
water (time, y, x) uint8 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
crs: EPSG:32650# Uncomment the following line to display WOfS for each time slice. This may take a few minutes
# wofl.water.plot(col="time", col_wrap=5);
Water observations summaries, based on odc-stats/plugins/wofs.py, consist of:
count_clear: a count of every time a pixel was observed (not obscured by terrain or clouds)count_wet: a count of every time a pixel was observed and wetfrequency: what fraction of time (wet/clear) was the pixel wet# Rename dimensions as required
wofl = wofl.rename({"x": "longitude", "y": "latitude"})
from odc.algo import safe_div, apply_numexpr, keep_good_only
wofl["bad"] = (wofl.water & 0b0111_1110) > 0
wofl["some"] = apply_numexpr("((water<<30)>>30)==0", wofl, name="some")
wofl["dry"] = wofl.water == 0
wofl["wet"] = wofl.water == 128
wofl = wofl.drop_vars("water")
for dv in wofl.data_vars.values():
dv.attrs.pop("nodata", None)
# Ignore warnings triggered by time slices without data at all
warnings.filterwarnings("ignore", message="divide by zero encountered in true_divide")
warnings.filterwarnings("ignore", message="invalid value encountered in true_divide")
wofl.wet.plot(col="time", col_wrap=5);
# Helper frunction from https://github.com/opendatacube/odc-stats/blob/develop/odc/stats/plugins/wofs.py
def reduce(xx: xr.Dataset) -> xr.Dataset:
nodata = -999
count_some = xx.some.sum(axis=0, dtype="int16")
count_wet = xx.wet.sum(axis=0, dtype="int16")
count_dry = xx.dry.sum(axis=0, dtype="int16")
count_clear = count_wet + count_dry
frequency = safe_div(count_wet, count_clear, dtype="float32")
count_wet.attrs["nodata"] = nodata
count_clear.attrs["nodata"] = nodata
is_ok = count_some > 0
count_wet = keep_good_only(count_wet, is_ok)
count_clear = keep_good_only(count_clear, is_ok)
return xr.Dataset(
dict(
count_wet=count_wet,
count_clear=count_clear,
frequency=frequency,
)
)
summary = reduce(wofl)
A count of every time a pixel was observed and wet.
summary.count_wet.plot(size=10)
plt.title(f"{aoi_label} – Wet observations counts for {time[0]} to {time[1]}");
A count of every time a pixel was observed (not obscured by terrain or clouds).
summary.count_clear.plot(size=10)
plt.title(f"{aoi_label} – Clear observations counts for {time[0]} to {time[1]}");
What fraction of the time was the pixel wet.
summary.frequency.plot(size=10)
plt.title(f"{aoi_label} - Wet frequency for {time[0]} to {time[1]}");
# Save frequencies to high-res image file
summary.frequency.plot(size=20)
plt.title(f"{aoi_label} - Wet frequency for {time[0]} to {time[1]}")
plt.savefig('/home/jovyan/tempe_wofs_2020.png', dpi=600);